home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / randist / exppow.c < prev    next >
Encoding:
C/C++ Source or Header  |  2001-05-14  |  3.5 KB  |  133 lines
  1. /* randist/exppow.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <math.h>
  22. #include <gsl/gsl_math.h>
  23. #include <gsl/gsl_sf_gamma.h>
  24. #include <gsl/gsl_rng.h>
  25. #include <gsl/gsl_randist.h>
  26.  
  27. /* The exponential power probability distribution is  
  28.  
  29.    p(x) dx = (1/(2 a Gamma(1+1/b))) * exp(-|x/a|^b) dx
  30.  
  31.    for -infty < x < infty. For b = 1 it reduces to the Laplace
  32.    distribution. 
  33.  
  34.    The exponential power distribution is related to the gamma
  35.    distribution by E = a * pow(G(1/b),1/b), where E is an exponential
  36.    power variate and G is a gamma variate.
  37.  
  38.    We use this relation for b < 1. For b >=1 we use rejection methods
  39.    based on the laplace and gaussian distributions which should be
  40.    faster.
  41.  
  42.    See P. R. Tadikamalla, "Random Sampling from the Exponential Power
  43.    Distribution", Journal of the American Statistical Association,
  44.    September 1980, Volume 75, Number 371, pages 683-686.
  45.    
  46. */
  47.  
  48. double
  49. gsl_ran_exppow (const gsl_rng * r, const double a, const double b)
  50. {
  51.   if (b < 1) 
  52.     {
  53.       double u = gsl_rng_uniform (r) ;
  54.       double v = gsl_ran_gamma (r, 1/b, 1.0) ;
  55.       double z = a * pow(v, 1/b) ;
  56.  
  57.       if (u > 0.5) 
  58.     {
  59.       return z ;
  60.     } 
  61.       else 
  62.     {
  63.       return -z ;
  64.     }
  65.     }
  66.   else if (b == 1) 
  67.     {
  68.       /* Laplace distribution */
  69.       return gsl_ran_laplace (r, a) ;
  70.     }
  71.   else if (b < 2) 
  72.     {
  73.       /* Use laplace distribution for rejection method */
  74.  
  75.       double x, y, h, ratio, u ;
  76.  
  77.       /* Scale factor chosen by upper bound on ratio at b = 2 */
  78.  
  79.       double s = 1.4489 ; 
  80.       do 
  81.     {
  82.       x = gsl_ran_laplace (r, a) ;
  83.       y = gsl_ran_laplace_pdf (x,a) ;
  84.       h = gsl_ran_exppow_pdf (x,a,b) ;
  85.       ratio = h/(s * y) ;
  86.       u = gsl_rng_uniform (r) ;
  87.     } 
  88.       while (u > ratio) ;
  89.       
  90.       return x ;
  91.     }
  92.   else if (b == 2)
  93.     {
  94.       /* Gaussian distribution */
  95.       return gsl_ran_gaussian (r, a/sqrt(2.0)) ;
  96.     }
  97.   else
  98.     {
  99.       /* Use gaussian for rejection method */
  100.  
  101.       double x, y, h, ratio, u ;
  102.       const double sigma = a / sqrt(2.0) ;
  103.  
  104.       /* Scale factor chosen by upper bound on ratio at b = infinity.
  105.      This could be improved by using a rational function
  106.      approximation to the bounding curve. */
  107.  
  108.       double s = 2.4091 ;  /* this is sqrt(pi) e / 2 */
  109.  
  110.       do 
  111.     {
  112.       x = gsl_ran_gaussian (r, sigma) ;
  113.       y = gsl_ran_gaussian_pdf (x, sigma) ;
  114.       h = gsl_ran_exppow_pdf (x, a, b) ;
  115.       ratio = h/(s * y) ;
  116.       u = gsl_rng_uniform (r) ;
  117.     } 
  118.       while (u > ratio) ;
  119.  
  120.       return x;
  121.     }
  122. }
  123.  
  124. double
  125. gsl_ran_exppow_pdf (const double x, const double a, const double b)
  126. {
  127.   double p ;
  128.   double lngamma = gsl_sf_lngamma (1+1/b) ;
  129.   p = (1/(2*a)) * exp(-pow(fabs(x/a),b) - lngamma);
  130.   return p;
  131. }
  132.  
  133.